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SUMMARY 


^Differential  equations  governing  the  transient  response  of  thermal  protec¬ 
tion  systems  to  a  hyperthermal  environment  are  presented.  These  equations  are 
expanded  into  finite-difference  equations  which  are  suitable  for  numerical  solu¬ 
tion.  The  equations  provide  for  three  layers  of  different  materials,  the  first 
two  of  which  may  have  moving  boundaries.  Concentrated  heat  sinks,  such  as 
metallic  structures,  may  be  located  at  the  back  surfaces  of  the  second  or  third 

layers  or  of  both  layers ^  \ 

-  }  U^-<r 

The  analysis  was  developed  primarily  for  charring  ablators  but  is  also 
applicable  to  impregnated  ceramic,  subliming,  and  heat-sink  thermal  protection 
systems.  The  principal  difficulty  encountered  in  numerical  analysis  of  charring 
ablators  is  the  extensive  computer  time  required  to  obtain  solutions.  Provision 
is  made  in  the  numerical  equations  to  introduce  options  which  reduce  computer 
time.  The  errors  resulting  from  these  options  under  various  conditions  are 
discussed.  Good  agreement  is  obtained  between  numerical  results  and  exact 
solutions \  v 


INTRODUCTION 


An  adequate  thermal  protection  system  may  constitute  20  to  30  percent  of 
the  total  reentry  weight  for  vehicles  which  must  enter  the  earth1 s  atmosphere 
at  supercircular  velocity.  Preliminary  studies  (ref.  1)  indicate  that  charring 
ablators  provide  the  most  efficient  thermal  protection  shield  for  most  of  the 
vehicle.  Other  materials  may  be  required  for  certain  vehicle  areas. 

Extensive  experimental  investigations  have  been  conducted  on  the  perform¬ 
ance  of  charring  ablators.  (See  refs.  2  to  6.)  Most  of  this  work  is  concerned 
with  the  overall  performance  of  the  material  in  a  given  environment.  The  theo¬ 
retical  analysis  of  charring  ablator  systems  is  very  complex  by  comparison  with 
the  analysis  of  other  heat-shield  systems,  and  many  approximations  and  assump¬ 
tions  must  be  made  to  obtain  solutions.  Several  procedures  for  the  numerical 
analysis  of  the  performance  of  these  materials  during  reentry  have  been  reported. 
These  analyses  differ  primarily  in  the  treatment  of  thermal  decomposition 


processes.  In  the  more  rigorous  approach  (refs.  7  and  8),  details  of  the  chem¬ 
ical  kinetics  of  the  reactions  through  the  depth  of  the  material  are  considered. 
In  the  more  simplified  approach  to  the  problem  (ref.  9),  it  is  assumed  that  the 
decomposition  processes  occur  in  a  single  plane  at  a  temperature  which  may  be 
fixed  or  a  function  of  the  rate  of  decomposition.  Since  it  has  been  shown  that 
the  thickness  of  the  decomposition  zone  is  very  small  compared  with  the  total 
thickness  of  the  material  (ref.  10),  this  simplifying  assumption  appears  to  be 
valid.  In  the  present  analysis,  the  simplified  approach  is  followed. 

The  most  difficult  problem  encountered  in  analysis  of  the  performance  of 
charring  ablators  is  the  formulation  of  a  quantitative  expression  for  the  rate 
of  char  removal.  It  appears  that  the  char  may  be  removed  chemically  (oxida¬ 
tion),  thermally  (sublimation),  mechanically  (spalling,  aerodynamic  shear),  or 
by  a  combination  of  these  methods,  depending  on  the  specific  material  involved. 
Results  of  an  investigation  of  oxidation  as  a  mechanism  in  the  performance  of 
charring  ablators  are  presented  in  reference  11.  A  theoretical  model  which  pre¬ 
dicts  spalling  of  the  char,  which  is  observed  in  some  tests,  is  presented  in 
reference  12.  It  is  also  anticipated  that  aerodynamic  shear  forces  have  a  pro¬ 
nounced  effect  on  the  removal  of  low-density  chars.  An  analytical  program  must 
have  considerable  flexibility  if  it  is  to  be  suitable  for  investigation  of  the 
various  mechanisms  of  char  removal  that  may  be  operative. 

In  this  paper,  equations  are  derived  for  calculating  the  thermal  response 
of  charring  ablators  to  reentry  heating  conditions.  Sufficient  options  are  pro¬ 
vided  so  that  not  only  may  all  methods  of  char  removal  be  considered  but  also, 
by  proper  manipulation,  the  response  of  impregnated  ceramics,  subliming  abla¬ 
tors,  and  heat- sink  materials  to  severe  heating  conditions  can  be  determined. 

The  equations  have  been  programed  for  solution  by  an  electronic  data  processing 
system,  and  numerical  examples  are  presented  to  demonstrate  the  accuracy  of 
various  approximations.  The  equations  presented  here  are  similar  to  those  in 
reference  9*  The  primary  difference  is  that  in  this  paper  the  finite-difference 
equations  are  derived  for  fixed  points  in  a  moving  coordinate  system  whereas  in 
reference  9  moving  points  in  a  fixed  coordinate  system  are  used.  The  present 
paper  also  provides  greater  flexibility  in  the  boundary  conditions  which  can  be 
imposed. 


SYMBOLS 


The  units  used  for  the  physical  quantities  defined  in  this  paper  are  given 
both  in  the  U.S.  Customary  Units  and  in  the  International  System  of  Units  (SI). 
Factors  relating  the  two  systems  are  given  in  reference  13. 

A  constant  in  mass  loss  rate  equation  corresponding  to  specific 

reaction  rate 

B  constant  in  exponential  of  mass  loss  rate  equation  corresponding  to 

activation  energy 
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b  thickness  of  layer  at  which  transfer  is  made  between  complete  equa¬ 

tions  and  boundary  equations  only  (see  appendix  C) 

C0  oxygen  concentration  by  weight  external  to  boundary  layer 

Ci+j,  Ci+j+m  heat  capacity  of  heat  sink  (Pcpx)heat  sink 

Cw  oxygen  concentration  by  weight  at  wall 

Ci,  C2  constants  of  integration  (eq.  (5^)) 


D 

F 

G 

g 

Hc 

he 

hw 

£hc 

Ahf 

Ahp 

i 

A 

K 

k 

kn,  n+1 


l 


specific  heat 

specific  heat  of  gaseous  products  of  pyrolysis 
parameter  defined  by  equation  (53)  j  j(%^p  +  mpcp) 

heat-generation  function 

constant  in  equation  (9d) 

initial  temperature  distribution 

heat  of  sublimation  of  char 

local  enthalpy  external  to  boundary  layer 

local  enthalpy  of  fluid  at  wall 

heat  of  combustion 

heat  capacity  of  coolant 

heat  of  pyrolysis 

number  of  stations  in  char,  including  one  at  each  surface 
number  of  stations  in  uncharred  material,  including  one  at  backface 
reaction- rate  constant  for  oxidation  of  char 
thermal  conductivity 

thermal  conductivity  evaluated  at  temperature  of  point  midway  between 
stations  n  and  n+1 

distance  between  stations 
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in 


m 


mn 


mr 


m, 


ox 
NLe 
Pw 
q 


q, 


aero 


qG, net 

% 


R 

S(0  -  T) 
T 


T- 


B 


T- 


n 


Ti 


number  of  stations  in  insulation,  including  one  at  back  surface 
mass  loss  rate,  acmc  +  ccpip 

rate  of  char  loss 

rate  of  loss  of  uncharred  material 

rate  at  which  oxygen  diffuses  to  surface 

Lewis  number 

total  pressure  at  wall 

heating  rate 

net  aerodynamic  heating  rate  at  the  surface  (eq.  (25b)) 
convective  heating  rate  to  a  cold  wall  with  no  mass  transfer 
hot-wall  convective  heating  rate  corrected  for  transpiration 
radiant  heating  rate 
effective  nose  radius 

unit  step  function:  S  =  0,  0  <  Tj  S  =  1,  0^T 

temperature 

temperature  to  which  back  surface  radiates 
temperature  at  finite-difference  stations 
temperature  of  pyrolysis 


^i+j>  -^i+j+m  temperature  at  which  cooling  system  is  activated 
T-]_  temperature  at  which  char  sublimes 

t  time 

V  velocity 
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AWf  rate  of  coolant  consumption 

x  thickness 

x  distance  from  initial  outer  surface  to  outer  surface  of  char  layer 

y  distance  from  initial  outer  surface  of  char  layer 

Z  distance  from  back  surface  of  shield 

a  absorptivity  of  char  surface 

ac  weighting  factor  for  transpiration  effectiveness  of  char  mass  loss 

Op  weighting  factor  for  transpiration  effectiveness  of  pyrolysis  products 

(3  either  0  or  1  depending  on  whether  transpiration  or  ablation  theory  is 

used 

e  emissivity 

£  transformed  coordinate  for  uncharred  material  (eq.  (23b)) 

T)  transpiration  coefficient 

0  temperature  in  differential  equations 

A  weight  of  char  removed  per  unit  weight  of  oxygen 

|  transformed  coordinate  for  char  layer  (eq.  (23a)) 

p  density 

a  Stefan-Boltzmann  constant 

Subscripts: 

o  initial  value 

i, j,m  number  of  stations  in  first,  second,  and  third  layers  (see  fig.  2) 
n  integer 

s  stagnation  point 

oo  free  stream 

CALC  calculated  by  numerical  method 

EXACT  exact  solution 
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Unprimed  symbols  refer  to  char  layer  unless  otherwise  specified;  single 
primes  refer  to  uncharred  material;  double  primes  refer  to  insulation. 


ANALYSIS 


The  thermal  protection  system  that  is  to  be  analyzed  is  shown  schematically 
in  figure  1.  Although  this  discussion  is  confined  to  a  charring  ablator  system, 

all  the  concepts  and  equations  apply  equally  as  well  to  any 
other  thermal  protection  system  composed  of  not  more  than 
three  primary  layers .  For  a  charring  ablator  system,  the 
outer  (heated)  layer  is  the  char,  the  center  layer  contains 
the  uncharred  material,  and  the  third  layer  consists  of 
insulation.  Heat  sinks  can  be  located  at  the  back  of  the 
second  or  third  layers  or  at  both  locations. 

Char  Layer 


_ ^i2iaj_Cha_r_  Surface^  _ 


Uncharred  Material 


The  outer  (char)  surface  is  subjected  to  aerodynamic 
heating.  The  char  layer  provides  both  insulation  and  a 
high-temperature  outer  surface  for  reradiation.  The  heat 
passing  through  this  layer  is  partially  absorbed  by  pyrol¬ 
ysis  at  the  interface  between  the  char  layer  and  the 
uncharred  material,  and  the  remaining  heat  is  conducted  into 
the  uncharred  material.  The  gases  generated  by  pyrolysis 
transpire  through  the  char  layer  and  are  injected  into  the 
boundary  layer.  The  gases  are  heated  as  they  pass  through 
the  char,  and  this  heat  removal  from  the  char  layer  reduces 
the  quantity  of  heat  conducted  to  the  pyrolysis  interface. 
When  these  gases  are  injected  into  the  boundary  layer,  the 
convective  heat  transfer  is  reduced.  This  reduction  in  con¬ 
vective  heating  is  the  same  effect  as  that  obtained  with 
simple  subliming  ablators.  In  addition  to  the  gases  produced  by  pyrolysis,  the 
carbonaceous  residue  remaining  at  the  interface  adds  to  the  thickness  of  the 
char  layer.  While  the  processes  of  pyrolysis,  transpiration,  and  injection  are 
underway,  char  removal  may  also  be  taking  place  as  a  result  of  thermal,  chemical, 
or  mechanical  processes.  Thus  the  total  char  thickness  may  increase  or  decrease 
depending  on  the  relative  rates  of  formation  and  removal  of  the  char.  The  var¬ 
ious  processes  discussed  are  related  quantitatively  in  the  following  sections. 


Figure  1.-  Schematic 
diagram  of  system 
employing  charring 
ablator. 


Differential  Equations 


It  is  assumed  that  thermal  properties  in  a  given  layer  of  material  are 
functions  only  of  temperature,  that  all  heat  flow  is  normal  to  the  surface,  and 
that  gases  transpiring  through  the  char  are  at  the  same  temperature  as  the  char. 
Then  the  governing  differential  equations  (from  ref.  9)  for  the  char  layer 
(x^y^x  +  x)  are  as  follows: 
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(1) 


Heat  conducted  Heat  absorbed  by  Heat  generated  Heat  stored 
transpiring  gases 


for 

the 

uncharred  layer  (x  +  x  ^  y 

^  xo 

+  xi)> 

|=  p'c’  Ml 

dy  \ 

Sy  j 

>  p  at 

and 

for 

a  layer  of  insulation  ^xD 

+  x^ 

=  y  =  xD  +  x^ 

— (k" 

Ml 

\  =  p"c"  Ml 

Sy  \ 

Sy  j 

1  p  a 

+  x. 


(2) 


(3) 


The  thicknesses  of  the  layers  to  which  the  first  two  of  these  equations  apply 
vary  with  time  in  a  manner  which  is  determined  by  the  boundary  conditions. 


Initial  Conditions 

The  initial  temperature  distribution  is  assumed  to  be  given  as  a  function 
of  position: 


e(y,o)  =  g(y)  (4) 

The  initial  mass-transfer  rates  must  also  be  specified.  It  should  be  noted  that 
these  values  can  be  other  than  zero  for  some  cases. 


Surface  Boundary  Conditions 

Two  conditions  must  be  specified  at  the  heated  surface.  One  must  determine 
either  the  rate  of  removal  of  material  at  the  surface  or  the  temperature  of  the 
surface;  the  other  is  provided  by  the  energy  balance. 

Surface  ablation.-  In  general,  the  relative  importance  of  the  mechanisms 
involved  in  char  removal  from  specific  materials  is  not  well  established  at  this 
time.  It  has  been  established,  however,  that  oxidation  of  the  char  surface  is 
one  important  mechanism.  Spalling  of  the  char  as  a  result  of  internal  pressure 
is  observed  in  some  cases.  Ablation  at  a  given  temperature  (that  is,  sublima¬ 
tion)  occurs  if  the  heating  rate  is  sufficiently  high.  Ablation  of  the  surface 
may  also  occur  as  a  result  of  aerodynamic  shear  stresses. 

To  provide  maximum  flexibility,  provision  is  made  for  the  following  mecha¬ 
nisms  of  surface  erosion: 

(l)  Ablation  at  a  given  temperature  which  may  be  a  function  of  ablation 
rate  (sublimation) 


(2)  Removal  of  char  at  a  rate  which  is  a  given  function  of  time  (spalling, 
aerodynamic  shear) 

(5)  Removal  of  char  at  such  a  rate  that  the  char  thickness  is  a  given  func¬ 
tion  of  time  (spalling,  aerodynamic  shear) 

(4)  Ablation  as  a  result  of  a  chemical  process  (oxidation) 

For  ablation  at  a  given  temperature,  two  cases  are  considered.  In  one  case, 
ablation  occurs  at  a  fixed  temperature.  In  the  other  case,  the  char  mass  loss 
rate  is  an  exponential  function  of  the  surface  temperature.  For  ablation  at  a 
fixed  temperature,  no  surface  erosion  occurs  if  the  calculated  surface  tempera¬ 
ture  is  less  than  the  specified  ablation  temperature.  If  the  calculated  surface 
temperature  is  higher  than  the  ablation  temperature,  ablation  occurs  at  a  rate 
sufficient  to  reduce  the  temperature  to  the  ablation  temperature;  that  is,  me 
is  equal  to  zero  for  Tq  <  Tq,  and  me  is  calculated  from  an  energy  balance  at 
the  surface  for  Tq  =  Tq.  In  the  second  case,  char  mass  loss  rate  and  surface 
temperature  are  related  as  follows: 


mc  =  Ae-B/Tl  (5) 

An  equation  of  this  form  has  some  physical  signif icance,  because  decomposition 
reactions  proceed  more  rapidly  at  higher  temperatures.  By  an  appropriate  selec¬ 
tion  of  A  and  B,  equation  (5)  yields  results  similar  to  those  obtained  by 
specifying  an  ablation  temperature. 

If  the  rate  of  char  removal  is  given  function  of  time 

mc  =  f(t)  (6) 

then  the  rate  of  char  removal  is  obtained  from  the  input  data  and  the  surface 
temperature  is  calculated  from  an  energy  balance.  Such  a  relation  might  be  used 
to  compare  calculated  and  experimental  results  when  a  more  basic  quantitative 
relation  for  the  experimental  rate  of  char  removal  is  not  available. 

If  the  char  thickness  is  a  given  function  of  time 

x  =  f(t)  (7) 

then  the  rate  of  char  removal  is  calculated  from  this  relation  together  with  the 
rate  of  char  formation  which  is  calculated  from  the  conditions  at  the  pyrolysis 
interface.  This  condition  can  be  used  when  it  is  desirable  to  perform  calcula¬ 
tions  for  applications  in  which  the  char  thickness  is  known  as  a  function  of 
time  even  though  mechanisms  of  char  removal  may  be  present  which  cannot  be 
expressed  quantitatively. 

It  has  been  shown  experimentally  that  oxidation  is  an  important  mechanism 
of  char  removal.  (See  ref.  11. )  For  a  half-order  reaction,  the  rate  of  oxida¬ 
tion  of  carbon  can  be  determined  from  the  following  equation  (ref.  lk): 
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mr  =  Ae 


-B/T1^PW 


The  pressure  at  the  wall  must  he  specified  for  subsonic  and  supersonic  flow. 
However,  in  hypersonic  flow,  the  pressure  can  be  related  to  the  stagnation 
heating  rate  and  enthalpy.  The  stagnation  pressure  in  hypersonic  flow  is 
approximately  (see  ref.  15): 


Further  (from  ref.  16), 


p,  =  11  p  v2 

s  12  Mco  00 


a  oc  ^°°  y3 
q  v°° 


he  «  vf 


Then,  pw  s  can  be  expressed  by  the  relation 

=  GR(£l 


where  the  constant  of  proportionality  is 


G  .  4X0.72  ft3-sec,2-aH!  ,  56.2  m;-sec2-atlJ 


The  pressure  at  the  wall  is  therefore  given  hy 


PW  =  G  b(Jl 

pw,s  \hes 


where 


depends  on  vehicle  attitude  and  body  location.  The  rate  at  which 


char  is  removed  by  oxygen  must  be  proportional  to  the  net  rate  at  which  oxygen 
diffuses  to  the  surface.  From  reference  17  this  rate  is 


.  _  HLe  «C,net,  ,  _ 

"°x  '  ie  -  (°e  C»J  " 


(10a) 


As  shown  in  a  subsequent  section,  q^  is  the  hot-wall  convective  heating 

rate  corrected  for  transpiration  (see  eq.  (13));  that  is, 
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X^net  ~  h 


^  ~  - 11  - p)  °-724  ^(acAc +  v^) 

° '  15(^C )  (ac“c  +  “P^p)2  "  ^(“c^c  +  apAp)^ 

•d 


(10b) 


By  eliminating  the  concentration  of  oxygen  at  the  wall  Cw  in  equa¬ 
tions  (8)  and  (10a),  the  rate  of  removal  of  char  by  oxidation  is  found  to  be 


c  2 


1  )  (he  -  U)r2Pw  +  fee  +  4k2 

21  0.6  I  ^0.6  W 

^^net^Le  U^^net  Le 


(10c) 


where  K  =  Ae 


-B/Ti 


The  equation  for  a  first-order  oxidation  reaction  is  obtained  similarly. 
The  resulting  equation  is 


KPwCe 


(lOd) 


1  +  Kpw(he  ~  hw) 

,  0.6 

^C,net^Le 

Surface  location.-  When  char  removal  occurs ,  the  char  surface  moves  with 
respect  to  a  coordinate  system  fixed  in  the  material.  The  distance  between  the 
surface  of  the  char  and  the  initial  surface  location  is  given  by 


The  thickness  of  the  char  at  any  time  is  equal  to  the  initial  char  thickness, 
plus  the  thickness  of  char  formed  by  pyrolysis,  less  the  thickness  of  char 
removed;  that  is, 


x  =  xQ  + 


p1  -  p 


Surface  energy  balance.-  The  heat  input  consists  of  convective  and  radiant 
heating.  This  energy  must  be  accommodated  at  the  surface  by  a  combination  of 
four  mechanisms: 
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(1)  Blocking  by  mass  transfer  into  the  boundary  layer 

(2)  Reradiation  or  reflection  from  the  surface 

(3)  Conduction  into  the  material 

(4)  Sublimation  of  the  char 

The  effect  of  mass  transfer  on  heat  transfer  has  been  studied  extensively. 
With  low-mass-transfer  rates  it  is  found  that  the  reduction  in  heat-transfer 
rate  is  directly  proportional  to  the  product  of  the  mass-transfer  rate  and  the 
enthalpy  difference  across  the  boundary  layer.  With  high-mass-transfer  rates, 
which  may  occur  when  a  large  fraction  of  the  heat  input  is  radiant,  the  linear 
approximation  is  no  longer  adequate  and  it  is  necessary  to  use  a  higher  order 
approximation.  A  second-degree  approximation  is  derived  in  appendix  A. 

The  surface  energy  balance,  expressed  in  a  form  in  which  either  approxima¬ 
tion  to  the  blocking  effectiveness  can  be  selected,  is  as  follows: 


Cold  wall  Hot  wall 

convective  correction 

heating  rate 


Aerodynamic  blocking 


- - — ^ — - - 

Net  convective  heating  rate 


J 


+ 


aqR 


Jl  -  s(e  -  TjJjic  Ah, 


Radiative  Combust ive 

heating  rate  heating  rate 


k 

aelel 


Re radiation 


Conduction  Heat  of  sublimation 

to  interior  of  char 


(13) 


If  transpiration  theory  (second-degree  approximation,  appendix  A)  is  used, 

(3  is  equal  to  zero.  For  linear  ablation  theory,  3  is  equal  to  1.  In  either 
case,  the  heat  absorbed  by  vaporization  of  the  char  Hc  and  the  heat  of  com¬ 
bustion  of  the  char  Ahc  are  considered  separately.  The  coefficients  ac  and 
ap  can  be  used  to  differentiate  between  the  blocking  effectiveness  of  the  gases 
produced  at  the  surface  and  at  the  pyrolysis  interface.  Evaluation  of  these 
coefficients  is  discussed  briefly  in  appendix  A. 

The  heat  transfer  to  the  outer  surface  is  assumed  to  be  a  given  function 
of  time  and  consists  of  the  cold- wall  convective  heating  rate  q^  and  the  radi¬ 
ant  heating  rate  q^>  incident  on  the  surface.  These  two  components  must  be 

specified  separately  because  mass  transfer  at  the  surface  blocks  part  of  the 
aerodynamic  heating  but,  in  general,  has  no  effect  on  radiant  heating.  Addi¬ 
tional  terms  can  easily  be  included  in  equation  (13)  to  account  for  other 
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phenomena  which  may  affect  the  heat  balance  at  the  char  surface.  For  example, 
reference  18  discusses  a  gas-phase  combustion  in  the  boundary  layer  involving 
the  gases  of  pyrolysis.  This  effect  has  not  been  clearly  identified  at  the 
Langley  Research  Center  and  is,  therefore,  not  included  in  the  equation.  How¬ 
ever,  phenomena  such  as  this  may  be  important  in  some  cases  and  their  existence 
should  certainly  be  considered. 

Equation  (13)  is  normally  used  in  this  analysis  as  the  boundary  condition 
on  the  temperature  at  the  outer  surface.  However,  when  0  is  equal  to  the  sub¬ 
limation  temperature  Tq,  the  specified  sublimation  temperature  provides  the 
boundary  condition  on  the  temperature  and  equation  (13)  is  used  to  calculate  the 
rate  of  ablation  mc. 


Pyrolysis- Interface  Boundary  Condition 

Energy  balance.-  The  heat  conducted  to  the  pyrolysis  interface  must  be 
either  absorbed  by  pyrolysis  reactions  or  conducted  into  the  uncharred  material; 
that  is,  at  y  =  x  +  x, 


-k 


^  -  mp  Ahp 


(1U) 


In  addition,  the  temperatures  in  the  char  and  in  the  uncharred  material  must  he 
equal  at  the  interface;  that  is,  at  y  =  x  +  x, 

0=0'  (15) 

Pyrolysis  rate.-  Two  approaches  are  available  for  calculating  the  rate  of 
pyrolysis.  In  the  first  approach,  it  is  assumed  that  pyrolysis  occurs  at  a 
given  temperature  Tj_.  If 


o_  ^  rp 
DX+X  1 


then 


mp  =  0 


If 


( l6a) 


(l6b) 


ex+x  “  Ti 

the  temperature  is  known  and  equation  (lb)  is  used  to  calculate  the  rate  of 
pyrolysis. 

In  an  alternate  approach,  it  is  assumed  that  the  rate  of  pyrolysis  is  a 
known  function  of  temperature,  for  example 
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mp  =  A'e 


-B'/05+x 


(IT) 


when 


^x+x  ^ 

In  this  case,  equations  (l4)  and  (17)  are  solved  for  both  temperature  and  pyrol¬ 
ysis  rate.  The  value  of  Ti  is  still  specified,  and  if  this  temperature  is 
reached,  the  pyrolysis  rate  is  determined  only  from  equation  (l4)  so  that  this 
temperature  is  not  exceeded. 


Pyrolysis- interface  location.-  As  pyrolysis  occurs,  the  interface  between 
the  char  layer  and  the  uncharred  material  moves  with  respect  to  a  fixed  coor¬ 
dinate.  Its  distance  from  the  initial  char  surface  location  is 

x  +  x  =  x0  +  f  — -E —  dt  (l8) 

J  0  P  ~  P 

The  instantaneous  thickness  of  the  uncharred  material  is 


x f 


xi  - 


(19) 


Boundary  Conditions  at  Back  Surface  of  Ablation  Material 

A  number  of  conditions  can  be  imposed  at  the  back  surface  of  the  ablation 
material,  depending  on  whether  additional  insulation  is  provided  or  some  pro¬ 
vision  is  made  for  temperature  control.  Whether  insulation  is  used  or  not,  the 
ablation  material  may  be  attached  to  a  thermally  thin  plate  which  functions  as 
a  concentrated  heat  sink. 


Three-layer  system.-  If  insulation  is  used,  the  temperature  of  the  ablation 
material  is  equal  to  the  temperature  of  the  insulation  at  their  interface. 


xo+xi 


=  0. 


xo+x 


! 

o 


From  an  energy  balance  at  the  back  surface  of  the  ablation  material, 


(20a) 


-k 


»  Be'  _  p 

sy  -  c±+j 


St 


(20b) 


Two-layer  system.-  If  no  insulating  layer  is  used,  the  back  surface  can 
be  assumed  to  be  perfectly  insulated,  cooled  at  a  given  temperature,  or  may 
exchange  radiation  with  a  sink  of  known  temperature  in  the  interior  of  the 
structure.  An  energy  balance  yields  the  following  equation: 
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(21) 


-k' 


Be' 

Sy 


Ci+j 


+  ael+j 


The  temperature  at  which  the  cooling  system  is  activated  is  The  choice 

of  conditions  is  accomplished  by  making  the  inapplicable  terms  equal  to  zero 
(that  is,  Ci+j  =  0  and/or  Ti+j  >  6'  and/or  e1+j  =  Oj. 


Boundary  Condition  at  Back  Surface 
of  Insulating  Material 


If  an  insulating  material  is  used  behind  the  ablating  material,  the 
boundary  condition  at  the  back  surface  (y  =  x0  +  +  Xq )  is  similar  to 

equation  (2l);  that  is. 


-k" 


Be_ 

Sy 


ci+j 


Be” 
+m  Bt 


+  s 


(e”  -  Ti+j+m)^f  Ahf  +  dei+j+m[(e")4  -  T. 


(22) 


Transformation  of  Coordinates 

The  equations  derived  in  the  preceding  discussion  are  similar  to  those  pre¬ 
sented  in  reference  9*  In  reference  9>  these  equations  are  expressed  in  finite- 
difference  form  and  solved  in  a  fixed  coordinate  system.  To  maintain  a  fixed 
number  of  stations  in  layers  of  varying  thicknesses,  it  is  necessary  to  change 
the  locations  of  the  stations  and  to  interpolate  to  determine  the  temperatures 
at  the  new  locations  after  each  step  in  the  calculation.  This  procedure  not 
only  increases  the  time  required  to  perform  the  computations,  but  also  intro¬ 
duces  a  small  error  in  each  step  of  the  calculation.  This  difficulty  can  be 
eliminated  by  transforming  the  equations  to  a  coordinate  system  in  which  the 
finite-difference  stations  remain  fixed,  and  the  coordinates  themselves  move 
to  accommodate  the  changes  in  the  locations  of  the  surfaces  of  the  different 
materials. 

The  y-coordinate  can  be  transformed  to  £-  and  ^-coordinates  in  the  char 
and  uncharred  layers,  respectively,  by  using  the  following  equations: 


(23a) 


y  -  xo 


i  = 


mr 


p' 


x' 


dt 


Ik 


P 


(23b) 


In  this  coordinate  system  the  outer  surface  remains  fixed  at  |  =  0.  The  inter¬ 
face  is  located  at  £  =  1  in  the  char  and  at  £  =  0  in  the  uncharred  material. 
The  back  surface  of  the  uncharred  material  is  located  at  £  =  1.  A  number  of 
advantages  result  from  the  use  of  this  double  transformation.  First,  the  char 
always  extends  from  |  =  0  to  |  =  1.  Therefore,  the  temperatures  tend  to  be 
more  nearly  steady  state  than  would  be  the  case  with  a  coordinate  system  fixed 
at  the  surface  only.  A  second  advantage  is  the  positive  location  of  the  pyrol¬ 
ysis  interface.  A  similar  transformation  would  also  be  very  beneficial  in 
locating  the  center  of  the  reaction  zone  when  the  pyrolysis  reactions  are  con¬ 
sidered  in  detail.  Because  the  reaction  zone  is  typically  very  thin,  a  very 
fine  finite-difference  network  is  required  to  analyze  it.  With  transformations 
similar  to  those  here,  the  center  of  the  reaction  zone  can  be  located,  and  the 
fine  network  can  be  restricted  to  this  region  rather  than  covering  the  entire 
range  of  possible  reaction-zone  locations. 


In  the  transformed  coordinate  system,  equations  (l)  and  (2)  are  as  follows 
(for  the  char  layer  and  uncharred  layer,  respectively): 


x2 


\  ,  1  So 

me  +  ril  C  [  t(  rn  J  0 

j  x  H 

mccp  mpcp  +  ^ Ip  1  _  p  *0)  CP 

de 

bt 


+  F  =  pcp  ^  (24a) 


1  d  /vi  de *\ 

(x‘)2^ 


i  Vcp 

x'  p'  -  p 


(i  _  A =  o'c'  Ml 
1  M  p  P  a 


(24b) 


The  boundary  conditions  are  as  follows: 
At  g  =  0, 

S  (e 


q  =  ae-,04  -  |  M  +  mc 
Haero  1  x  C 


where 


%.ero  =  a(lR  +  let1  "  ^  “  ^l0*724  ^(ac&c  +  “p^p) 


Tx)  (hc  +  Ahc)  -  £bc 


he 


h 


(25a) 


0  15(^)  t0^  +  apAp)~  ~  PT1(ctc&c  +  aP%>)tjj  (25b) 


at  6-1,  6=0, 


0=9' 


and 


_k  d©  _  -a,  ah  _  ki  Ml 
mP  ^P  X' 


(26a) 


(26b) 
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and  at  £  =  1  for  only  two  layers,  the  condition  at  the  back  surface 
(y  =  xG  +  *q)  is 


When  three  layers  are  used,  the  conditions  at  the  hack  of  the  second  layer  are 

0*  =0"  (28a) 


and 


k*  89 *  p  89  '  n  it  89 
x’  ^  8t  8y 


(28b) 


If  three  layers  are  used,  the  condition  at  the  back  surface  (y  =  xQ  +  Xq  +  x£) 
is 


■'i+j+m  ^ 


i+j+m) 


-i+j+m 


(e")4  - 


Finite-Difference  Equations 

The  locations  of  the  finite-difference  stations  are  shown  in  figure  2. 


Figure  2.-  Location  of  finite- 
difference  stations . 
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The  methods  used  to  change  the  differential  equations  to  finite-difference  form 
are  given  in  appendix  B.  A  summary  of  the  finite-difference  equations  obtained 
by  these  methods  is  presented  here. 

The  equation  at  the  first  station  (eq.  (B25))  is 

>%2  + Wo  -  -  a<=[s(ti  -  hXHc + *»<=)  -  m 

/yrpn 

+  i(-llTi  +  i8t2  -  9T3  +  2T4)(mccp  +  mpcp)  +  |  F  =  pcp  |  —  (30) 

The  equations  for  the  first  layer  (eqs.  (B9)  and  (BIO))  are 


■^n-1  “  ^n 


-  kn,  n+1  ~  +  |K+1  -  3Tn  -  ^n-1  -  Tn+g) 


n-1,  n  i 


mccp  "*■  ™p^p 


n  -  1/  mpP  .  \ 

+  "  mcjCp 

and 


n-l,n  Tn~Y  ^  -  kn, n+1  —  ' ^  +  |K+1  +  3Tn  -  6ln.x  +  Tn_2) 


+  ZF  =  Zpc 


AT 

P  A t 


n 


(2  Is  n  %  i  -  2)  (31) 


Tn  “  Tn+1  ,  1, 


mccp  +  mpcp 


,n-lfV  •  L 
+  iTTlF7^  "  m°  p 


AT 


n 


2  k 


+  ZF  =  pcpZ  At 

The  equation  at  the  pyrolysis  interface  (eq.  (B30))  is 
,  Ti+1  -  Ti  ,  Ti  "  Ti-l\  .  „  .  A  ftr  .  CPP 


(n  =  i  -  1)  (32) 


M+l  Z 


1  ki-l,i  i 


+  ZF  +  mp/fCp  +_ 


P'  -  P) 


i(llTj.  -  l»u 


+  9*1-2  -  Zh.}) 


Cpp' 


P1  -  p[6 


+  l8ri+1  -  9Ti+2  +  2T1+5j 


-  2Ah- 


-  ^pcpZ  +  p  cpZ 


/At 


(33) 
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The  equations  for  the  second  layer  (eqs.  (B15)  and  (Bl6))  are 


k 


Tn-1  -  Tn  Tn-Tn+1 


h-l,n  1 


-  k 


n, n+1  1  1 


and 


•  t  n  T  Arn 

+  |(6Tn+1  -  5Tn  -  2Tn-l  -  T^g)1  +■  j-—  =  l  'P'Cp  ^ 


(i+l<n<i  +  j-2)  (3M 


k. 


1  ^n-1  "  -^n  v  t  ^n  ~  ^n+1 


n-l,n  1 • 


n, n+1  1 


+  K^n+1  +  3Tn  -  6Tn_x  +  Tn_2) 


i  +  j  -  n  “pP'cp  _  .  ATn 

j  p '  _  p  "  p  CP  At 


(n  -  i  +  J  -  1)  (35) 


The  equations  at  the  back  of  the  second  layer  are  for  two  layers 

(eg.  (B33)) 


ki+j-l,i+j 


-  t.uj  _  0c1+j(Ti+J  -  4)  -  s(Ti+j  ■  Ti+j)*f 


,  1 1'  c 

P  cp  -  +  Ci+jl-^— 


(3  6) 


and,  for  three  layers  (eq.  (B36)) 


k.,- 


i+j-1, i+j  1 1 


Ti+j-l  "  Ti+j  v"  Ti+j  Ti+j+l 

- - -  -K.4  .4  , 


^i+j^i+j+l  21 


I  I  I  l  ,  tl  It  l  n 

=  p  cn  tr  +  p  C-D  ^7  +  °j 


AT.  . 

\  1+J 


'P  ~2  '  K  TP  2  1+J/  At 

The  equation  for  the  third  layer  (eq.  (B19)) 


(37) 


,  11  Tn-1  “  Tn  v"  ^n  "  ^n+1  _  n  »i,  n  ^n 
kn-l,n  - f" - n,  n+1  p  p  At 


(38) 
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The  equation  at  the  back  surface  of  the  third  layer  (eq.  (B39))  is 


II 

Ir 

i+j+m-1,  i+j+m 


•^i+j+m-l  ~  -^i+j+m 


aei+j+m  (■'■i+j+m  “  ^b)  +  ^^i+j+m  -  ^i+j+m)^f  ^f 


,  f  ii  it  l"  ,  n  \^i+j+m 

P  CP  2  Cl+J+m)~At 


In  some  cases  it  is  desirable  to  use  the  so-called  thin- layer  options, 
presented  in  appendix  C,  because  of  the  considerable  saving  of  computer  time 
involved  in  their  use.  When  these  options  are  used,  the  layer  in  question  is 
divided  into  two  half-elements  and  the  interior  stations  in  the  layer  are  no 
longer  considered.  The  first  layer  is  considered  "thin"  when  x  is  less  than 
some  prescribed  thickness  b,  and  the  second  layer  is  considered  "thin"  when 
xT  is  less  than  some  prescribed  thickness  b*.  The  finite-difference  equa¬ 
tions  used  with  the  thin-layer  options  are  as  follows. 

At  the  front  surface,  method  1  gives  the  following  equation  when  x  ^  b 
(eq.  (Cl)) 


Ti  -  Tj_ 


1,  i  3T 


+  Wo  -  a€ 


1*1  -  mc  s(tx  -  T^Hc  +  £hc)  -  Ahc 


Tf  -  Tl 


(®cCp  +  ipCp^  +  |  F  =  pcp  | 


X  AT1 


The  equations  at  the  interface  are,  for  x  ^  b  and  x'  >  b'  (eq.  (C2)) 


.  Ti+1  -  Ti  .  pcP 

i+1  p  ai,l  +  xF  +  r-ip  ,l  Cp  +  pl  _  - 


(T1  -  Tl) 


ctp*  p  -V 

TUt[|(-11T1  +  l8Ti+l  -  5Ti+2  +  ^i+j)  -  '  (pcPx  +  p'ci!'); 


for  x  ^  b  and  x'  ^  b'  (eq.  (C3)) 


i  Ti+1  -  T^  T^  _  Tp  /_  Pcp  \  .  v 

—p —  ai,i -V +  + ■*  (<*  +  ^)(Ti  -  Ti) 


I  I 

2h-(T1+J  -  Tl)  -  2fihpJ  »  (pcp 


pc  +  p'ci^ 
P  P  /At 
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and  for  x  >  b  and  xT  Is  b'  (eq.  (C4)) 


»  Ti+j  "  Ti 


— fr 


Ti  -  Ti_! 

4-1,1  j 


^  +  ^r)|(11Ti  “  +  9Tl-2  "  ^-3) 


P  -  P 


(Ti+j  -  Ti)  -  2£hpj  =  (pcpl  +  p'c^x'j^i 


The  equation  at  the  back  surface  of  the  second  layer  for  two  layers 
(eq.  (C5))  is 


t  Tj_  -  T-j  +  j 


%i+J  x' 


«i+j(Ti+3  -  te)  -  s(Ti+j  -  Ti+j)^f 

=  (rj '  cp  +  ci+o)' 


and  for  three  layers  (eq.  (c6))  is 


'  Tj  -  T1+j  „ 

Ki^i+j  xt  “  Ki+j,i+j+l 


Ti+j  -  Ti+j-l 


I  „  t! 

t  •  X  t  t!  11  7 

,  1  r»  1  _  _l_  /s  r* 


=  PcPT+pcPT  +  Ci+j 


Calculation  of  Temperatures  at  Fixed  Points 

It  is  possible  to  calculate  the  temperature  at  any  number  of  fixed  points 
Zn  where  Z  is  the  distance  from  the  back  surface.  (See  fig.  2.)  The  tem¬ 
perature  at  is  calculated  as  follows  (by  using  linear  interpolation): 


When 


find  station  N  such  that 


0  <  Zn  ^  x" 


(i  +  j  +  m  -  N  )l"  ^Zn^(i  +  j+  m-  N-  l)z" 


Then 


T(zn)  =  TN  +  (TN  -  TN+1) 


Zn  -  (i  +  j  +  m  -  N )l" 
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When 


*"  ^  Zn  <  *'  +  x” 

find  N  such  that 

x"  +  (i  +  j  -  N)l '  ^  Zn  ^  x"  +  (i  +  j  -  N  -  1)2' 

Then 

Zn  -  x"  -  (i  +  j  -  N)2  ' 


T(zn)  =  %  +  (%  -  Tjj+i)- 


2 


(4  7) 


When 


x'  +  x"  ^  Zn  <  x  +  x'  +  x" 

find  N  such  that 

x'  +  x"  +  (i  -  W)2  ^  Zn  ^  x'  +  x"  +  (i  -  N  -  1)2 

Then 

/  \  /m  m  \zn  “  x'  -  x"  -  (i  -  W)2 

T(zn)  =  %  +  (TN  -  Tn+1) - - -  (48) 

If  the  outside  surface  moves  past  the  largest  Zn  because  of  ablation  of  the 
surface,  then  the  temperature  at  this  point  is  meaningless  and  is,  therefore, 
not  calculated. 


RESULTS  AND  DISCUSSION 


Comparisons  have  been  made  of  the  numerical  results  obtained  when  the 
various  options  presented  in  this  report  are  used.  Comparisons  with  exact  solu¬ 
tions  have  also  been  made.  The  results  of  these  studies  are  presented  and  dis¬ 
cussed  in  the  following  sections. 


Comparison  With  Exact  Solutions 

The  accuracy  of  the  numerical  analysis  is  determined  for  two  cases  for 
which  exact  solutions  are  available.  One  case  provides  an  exact  solution  for 
a  homogeneous  nonablating  heat-sink  material  subjected  to  a  suddenly  applied 
constant  heating  rate.  The  other  case  is  for  quasi- steady- state  ablation.  In 
both  cases,  constant  material  property  values  and  environmental  values  are  used 
and  reradiation  from  the  surface  is  not  considered. 
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Heat-sink  case.-  The  exact  solution  for  the  heat- sink  case  is  obtained 
from  equation  (All)  of  reference  19-  When  written  with  the  symbols  used  in 
this  paper,  the  equation,  is 


(49) 


The  values  of  q  and  material  properties  used  to  obtain  a  solution  to  equa¬ 
tion  (49)  are  given  in  table  I.  These  values  do  not  represent  a  real  material 
but  were  chosen  to  simplify  the  calculation. 


TABLE  I.-  VALUES  USED  FOR  HEAT-SINK  EXACT  SOLUTION 


Quantity 

Value  in 

U.S.  customary- 
units 

Value  in 

SI  units 

1 . 

100  Btu/ft2-sec 

1135  kW/m2 

t . 

0.1  ft 

3.048  cm 

k . 

0.01  Btu/ft-sec-°R 

62.4  w/m-°K 

p . 

10  lb/ft^ 

160  kg/m^ 

cp . 

0.1  Btu/lb-°R 

0.4184  kj/kg-°K 

To . 

530°  R 

294°  K 

When  a  comparison  is  made  between  surface  temperature  histories  calculated 
from  this  equation  and  calculations  made  with  the  finite-difference  equations, 

y 

the  following  results  are  obtained.  At  the  heated  surface  — - — -  =  0,  the 

x  +  x 

error  is  about  -0.4°  R  (-0.2°  K)  out  of  about  1863°  p  (1035°  K)  when  T0  is 
530°  R  (29^°  K). 

Figure  3  shows  a  comparison  between  the  finite-difference  solutions  and 
the  exact  solution  for  the  first  few  time  steps  of  the  calculations.  The  cal¬ 
culations  shown  in  figure  3(a)  use  a  time  interval  of  1/409 6  second.  It  can 
be  seen  that  after  about  1 6  time  steps  (1/256  second  each),  the  agreement  is 
adequate. 

The  results  obtained  when  the  time  interval  is  cut  in  half  are  shown  in 
figure  3(h).  It  can  be  seen  that  the  error  in  surface  temperature  is  reduced 
considerably  for  the  same  heating  period. 
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(a)  At  =  1/4096  second. 


(b)  At  5=  1/8192  second. 


Figure  3*-  Comparison  of  surface  temperatures  for 
initial  time  steps. 


Quasi-steady-state  ablation  case.-  A  quasi-steady-state  case  is  one  in 
which  the  front  surface  and  pyrolysis  interface  recede  at  the  same  ratej  thus, 
the  char  thickness  is  constant.  If  conditions  are  also  chosen  such  that  there 
is  no  heat  transfer  into  the  uncharred  layer,  an  exact  solution  can  be  obtained. 


The  transformed  heat  balance  equation  for  the  first  layer  is  (eq.  (2^a)) 


JL_  A 
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With  quasi-steady-state  ablation,  x,  ip,  and  mc  are  constant.  Further, 


- - mc  =  0 

P  -  P 


Therefore  equation  (2ka)  becomes  (with  F  =  0): 


1  d  /.  d0  \  ,  1  d6  /  •  „  .  •  -  \  _  n 
x2  d|  \  d£  /  x  d|  V  c  P  Py 


Then  with  constant  k,  Cp  and  Cp,  equation  (51)  can  be  written 


H  +  d  —  =  0 

d£ 


where 


D  »  r-  + 


k  r^P 


“p5p) 


The  solution  of  equation  (52)  is  (0  being  replaced  by  T) 


T  =  Cp  +  C^e 


and  Cp  and  Cp  can  be  found  from  the  boundary  conditions: 

T  =  Tc 


U  -  0) 


T  =  Tt 


(I  =  1) 


The  final  solution  is 


TP  "  Tce~P  +  (TC  -  Tp)« 


1  -  e' 


To  determine  the  mass  loss  rates  mc  and  rbp,  and  the  char  thickness  x 
associated  with  this  equation,  the  following  equations  are  used: 


2k 


(56) 
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In  the  last  equation,  it  is  assumed  that  — 


^=0 


is  negligibly  small;  that  is, 


there  is  no  heat  flow  into  the  second  layer.  By  solving  equations  (55)  and 
(56)  for  mc>  Ap,  and  x,  the  following  equations  are  obtained: 


The  values  of  qaero  and  material  properties  used  to  obtain  a  solution 
to  equation  (55)  (by  using  eqs.  (55)  and  (57))  are  given  in  table  II.  These 
values  do  not  represent  a  real  material  but  were  chosen  to  simplify  the 
calculation. 


TABLE  II.-  VALUES  USED  FOR  QUASI  -  STEADY  -  STATE  EXACT  SOLUTION 


Quantity 

Value  in 

U.S.  Customary  Units 

Value  in 

SI  Units 

l . 

69.3  Btu/ft2-sec 

785  kW/m2 

k . 

1.0  X  10"^  Btu/ft-sec-°R 

0.62k  W/m-°K 

P . 

20  lb/ft3 

320  kg/m3 

p’ . 

bO  lb/ft^ 

64 0  kg/iu3 

fP . 

0.5  Btu/lb-°R 

2.09  kJ/kg-°K 

fi>  •  * . 

0.5  Btu/lb -°R 

2.09  kj/kg-°K 

Tl . 

4000°  R 

2222  °K 

Ti . 

1000°  R 

556  °K 

Ahp . 

1000  Btu/lb 

2.324  Mj/kg 

He . . 

1000  Btu/lb 

2.324  Mj/kg 
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(a)  Mass  loss  rate. 


(b)  Temperature. 

Figure  4.-  Error  in  quasi- steady-state  solu¬ 
tion  (x  =  0.01  ft  -  0.3048  cm). 


The  results  of  a  comparison 
between  the  quasi- steady- state  exact 
solution  and  the  numerical  analysis 
are  shown  in  figure  4.  For  the  case 

considered,  the  ratio  — - -  =  1 

P  P 

so  that  mc  =  mp.  Figure  4(a)  shows 

the  ratio  of  the  mass  loss  rate 
obtained  from  the  numerical  solution 
to  the  value  obtained  from  the  exact 
solution  plotted  as  a  function  of  the 
distance  l  between  stations.  In 
the  first  material,  the  error  is  less 
than  0*7  percent  even  with  large  l. 

The  ratio  of  the  temperatures 
obtained  from  the  numerical  solution 
to  those  obtained  from  the  exact  solu¬ 
tion  plotted  as  a  function  of  the  dis¬ 
tance  between  stations  in  the  first 
material  is  shown  in  figure  4(b)  for 
three  locations  through  the  thickness. 
The  error  is  greatest  at  the  center 
location,  as  would  be  expected,  since 
the  temperature  is  fixed  at  the  sur¬ 
face  and  interface .  Note  that  the 
error  is  less  than  0.2  percent  even 
for  large  l . 


Thin- Layer  Options 

The  time  interval  for  which  the  numerical  solution  is  stable  is  propor¬ 
tional  to  the  square  of  the  distance  between  finite-difference  stations.  (See 
appendix  C.)  In  general,  with  the  present  formulation  for  the  first  deriva¬ 
tives,  the  distance  between  stations  is  less  than  or  equal  to  one-fourth  the 
thickness  of  the  layers.  In  many  cases,  the  initial  value  of  x,  or  the  final 
value  of  x1,  is  very  small  and  excessive  computer  time  is  required.  To  reduce 
computer  time  in  such  cases,  equations  which  eliminate  the  interior  stations 
when  the  thickness  of  a  layer  is  less  than  a  specified  value  are  provided. 

These  equations  are  derived  in  appendix  C.  Specifically,  when  x  ^  b  (referred 
to  as  the  b-option),  the  interior  stations  in  the  first  layer  are  not  considered 
and  the  entire  layer  becomes  two  half  elements  of  thickness  x/2.  If  the  first 
layer  reaches  a  reasonable  thickness  due  to  pyrolysis  at  the  interface,  the 
thin-layer  option  is  no  longer  appropriate  and  the  program  reverts  to  the  stand¬ 
ard  equations  when  this  value  of  b  is  exceeded.  A  b-value,  which  is  one 
order  of  magnitude  larger  than  the  initial  value  of  x,  generally  provides 
acceptable  accuracy  and  significantly  reduces  computer  time.  Similarly,  when 
xr  ^  bT  (bf-option),  the  interior  stations  of  the  second  layer  are  discarded. 

A  b!-value,  one  order  of  magnitude  less  than  the  initial  value  of  x',  is  gen¬ 
erally  acceptable. 


2  6 


Results  obtained  with  the  thin-layer  options  are  compared  with  the  standard 
solution  in  figures  5  and  6  for  a  typical  charring  ablator  subjected  to  a  con¬ 
stant  heating  rate.  The  values  used  in  the  calculations  are  given  m  table  III. 
These  values  are  representative  of  a  low-density  phenolic-nylon  charring 
ablator. 


(a)  Char  layer. 


(h)  Uncharred  layer. 

Figure  5.-  Error  in  calculated  thickness  when  thin- char- 
layer  option  is  used. 
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TABLE  III.-  VALUES  USED  FOR  THIN- LAYER  OPTION  CALCULATIONS 


Quantity 

Value  in 

U.S.  Customary  Units 

Value  in  Si  Units 

4c . 

100  Btu/ft^-sec 

1135  kW/m2 

k  for  - 

500°  R  (278°  K) . 

5  X  10" 6  Btu/ft2-sec-°R 

0.0312  W/m-°K 

1000°  R  (556°  K)  . 

2  x  10"5  Btu/ft2-sec-°R 

0.1248  W/m-°K 

1500°  R  (834°  K)  . 

5  X  10“5  Btu/ft2-sec-°R 

0.312  W/m-°K 

2000°  R  (1110°  K)  . 

1.0  X  10-1*-  Btu/ft2-sec-°R 

0.624  W/m-°K 

7000°  R  (389O0  K)  . 

1.0  x  10~3  Btu/ft2-sec-°R 

6.24  W/m-°K 

CP . 

0.5  Btu/lb-°R 

2.09  kj/kg-°K 

5p . 

0.5  Btu/lb-°R 

2.09  kj/kg-°K 

Cp . 

0.32  Btu/lb-°R 

1.34  kj/kg-°K 

k' . 

3.06  X  10“5  Btu/ft-sec-°R 

0.191  w/m-°K 

P . 

15  lb/ft? 

240  kg/m3 

P' . 

36  lb/ft3 

576  kg/m3 

Ahp . 

250  Btu/lb 

O.581  MI /kg 

A . 

6.73  X  108  lb/ft2-sec-atm1/2 

3.29  X  109  kg/m2-sec-atm1/2 
2.216  X  XT  °K 

B . 

3.9877  X'lcA  °R 

ce . 

0.23 

X . 

0.75 

Pw,  atm  . 

1.3 

7.74  X  108  kg/m2-sec-atm'L/2 

A* . 

I.585  x  108  lb/ft2-sec- atm1/ 2 

B* . 

2.321  X  10^  °R 

1.29  x  left  °K 
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Figure  5(a)  shows  the  ratio  of 
char- layer  thickness  using  the  thin- 
char- layer  option  to  the  char-layer 
thickness  using  the  standard  equations. 

As  can  be  seen  from  the  figure,  when 
x  <  b  throughout  the  calculation,  there 
is  an  appreciable  error.  However,  for 
both  calculations,  the  char  thickness  is 
monotonically  increasing  and  if  a  b-value 
is  used  such  that  x  becomes  greater 
than  b  during  the  calculation,  the  char 
thickness  immediately  begins  to  approach 
the  thickness  obtained  when  the  option  is 
not  used.  Thus,  since  most  of  the  com¬ 
puter  time  necessary  for  a  calculation 
is  usually  consumed  in  the  beginning  of 
a  calculation  when  the  char  thickness  is 
small,  the  use  of  an  intermediate  b-value 
usually  produces  a  considerable  saving  of 
computer  time  without  any  appreciable 
error. 

Figure  5(b)  shows  the  ratio  of  the 
thickness  of  the  uncharred  layer  using 
the  thin- charred- layer  option  to  the 
uncharred  layer  thickness  using  the 
standard  equations.  The  error  which 
occurs  when  the  b-option  is  used  is  sig¬ 
nificant.  However,  as  with  the  char  layer 
reduces  this  error  considerably. 


(a)  Char  layer. 


(b)  Uncharred  layer. 

Figure  6.-  Error  in  calculated  thick¬ 
ness  vhen  thin-uncharred- layer 
option  is  used. 


the  use  of  an  intermediate  b-value 


Figure  6  shows  the  results  obtained  by  using  the  bT -option.  The  value  of 
b*  used,  0.02  inch  (0.05  cm),  is  such  that  x1  becomes  less  than  b1  at  about 
130  seconds.  The  error  shown  in  the  figure  can  be  reduced,  if  necessary,  by 
using  a  smaller  b* -value. 


CONCLUDING  REMARKS 


Differential  equations  governing  the  thermal  performance  of  charring  abla¬ 
tors  are  derived.  These  equations  are  expanded  into  finite-difference  form  and 
have  been  programed  for  numerical  solution  on  a  digital  computer.  Numerical 
results  compare  favorably  with  available  exact  solutions. 

The  required  computer  time  presents  a  major  problem  in  the  numerical  solu¬ 
tion  of  these  equations.  The  governing  equations  have  been  examined  and  a  num¬ 
ber  of  approximations  which  reduce  the  computer  time  are  introduced.  The  con¬ 
ditions  under  which  these  simplifications  should  be  used  and  the  error  involved 
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in  their  use  are  discussed.  The  computer  program  based  on  the  equations  pre¬ 
sented  herein  has  been  found  to  provide  a  practical  basis  for  heat-shield 
design  studies. 


Langley  Research  Center, 

National  Aeronautics  and  Space  Administration, 

Langley  Station,  Hampton,  Va.,  April  1,  1965* 
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APPENDIX  A 


EFFECTS  OF  MASS  INJECTION  ON  AERODYNAMIC  HEATING 


The  aerodynamic  heat  input  consists  of  both  convective  heating  and  radia¬ 
tive  heating.  The  material  that  is  liberated  by  erosion  of  the  outer  surface 
and  pyrolysis  at  the  interface  enters  the  boundary  layer  and  causes  a  signifi¬ 
cant  reduction  in  convective  heating.  It  has  been  shown  (refs.  20,  21,  and  22) 
that,  for  moderate  mass-transfer  rates,  the  effectiveness  of  mass  transfer  in 
blocking  convective  heating  is  approximately  a  linear  function  of  the  enthalpy 
in  the  flow  outside  the  boundary  layer.  It  is  further  shown  in  references  22 
and  23,  however,  that  this  approximation  may  lead  to  serious  error  if  the  sur¬ 
face  is  exposed  to  a  net  radiant  heat  input  such  as  that  which  a  typical  vehi¬ 
cle  would  experience  from  reentry  at  parabolic  or  hyperbolic  velocities. 

Blocking  effectiveness  as  a  function  of  a  mass-transfer-rate  parameter 
mhe^qQ  ^ where  m  =  acmc  +  OpHip)  is  shown  in  figure  7*  The  exact  solution  was 

obtained  from  the  boundary-layer  solutions  for  air-to-air  injection  (ref.  21). 
It  can  be  seen  from  figure  7  that  the  usual  linear  approximation  for  a  laminar 
boundary  layer  is  not  valid  at  higher  values  of  the  mass-transfer-rate  param¬ 
eter.  The  effect  of  radiant  heating  is  to  increase  m  and  this  increase 
results  in  operation  at  higher  values  of  mhg/q^.  The  second-degree  approxi¬ 
mation  was  developed  by  fitting  a  curve  at  values  of  mheyqQ  of  0,  1.0, 
and  2.5. 


Very  little  experimental  data  is  available  at  high  mass-transfer  rates, 
and  in  conservative  calculations  it  is  desirable  to  specify  a  minimum  value  of 


^Ghnet 

- r  greater  than  0. 

M 

he  ) 

This  approach  is  certainly 
in  order  if  boundary- layer 
separation  is  found  to 
occur  at  high  mass- 
injection  rates.  In  this 
analysis,  a  minimum  value 
of  the  q- ratio  (0.04-)  was 
specified  where  the  second- 
degree  approximation  departs 
from  the  theoretical  curve. 

The  equation  for  the 
second-degree  approximation 
is  as  follows: 


Figure  Blocking  effectiveness  for  a  laminar  bound¬ 
ary  layer  -with  air-to-air  injection. 
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APPENDIX  A 


The  constants  ac  and  Op  are  used  to  correct  equation  (Al)  for  the  differ¬ 
ence  between  the  molecular  weight  of  the  boundary- layer  gas  and  the  molecular 
weight  of  the  injected  gas.  The  constant  ac  must  also  be  corrected  for  that 
part  of  the  char  that  is  removed  mechanically,  that  is,  for  the  char  that  is 
removed  but  not  sublimed.  A  similar  approach  might  be  used  to  correct  for 
turbulent  flow.  The  effects  of  molecular  weight  and  turbulent  flow  are  dis¬ 
cussed  in  reference  2k.  An  equation  of  the  form  of  equation  (Al)  is  used  to 
correlate  experimental  data  in  reference  18  by  using  the  correction  of  refer¬ 
ence  2^4-  for  molecular  weight  effects.  The  equations  have  been  formulated  to 
provide  the  option  of  using  either  the  linear  function  of  enthalpy  (ablation 
theory)  or  a  second-degree  approximation  to  the  actual  blocking  effectiveness 
(transpiration  theory) . 
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APPENDIX  B 


FINITE-DIFFERENCE  EQUATIONS 


The  system  of  partial  differential  equations  to  be  solved  is  extremely 
complicated  (nonlinear  with  variable  coefficients)  and  solutions  in  analytic 
form  can  hardly  be  expected.  The  systems  can  be  put  into  a  form  suitable  for 
numerical  calculation  by  deriving  the  equations  in  finite-difference  form  and 
then  solving  them  by  an  iteration  procedure  on  a  high-speed  digital  computer. 


For  this  type  of  procedure  it  is  convenient  to  space  the  finite-difference 
stations  at  equal  intervals  throughout  a  given  layer,  as  shown  in  figure  2. 

The  first  (char)  layer  contains  i  stations  including  one  at  the  front  and  one 
at  the  back  surface  of  the  first  layer.  The  second  layer  contains  j  stations 
including  one  at  the  back  surface  of  the  second  layer.  The  third  layer  (if 
desired)  contains  m  stations  including  one  at  the  back  surface. 


The  distances  between  stations  are  therefore 


-\ 


(Bl) 


for  the  three  layers  whose  thicknesses  are  x,  x’,  and  x",  respectively.  In 
this  appendix,  the  finite-difference  form  of  the  general  equations  (eqs.  (24a), 
(24b),  and  (3))  are  presented  first  and  then  the  finite-difference  form  of  the 
boundary  equations  (eqs.  (25a),  (26b),  (27),  and  (29))  are  given. 


Interior  Stations 


Char  layer  (2  ^ 


- 1-4 12)  - 

X2  5|  \  5|  J 


n  %  i  -  l).-  The  finite-difference  form  of  equation 


1  59 
x  5| 


mccp  +  mpCp  +  | 


“pP 


He  C, 


+  F  =  pcp  — 


5e 

5t 


( 24a) 
(B2) 


is  obtained  as  follows:  The  second  derivative  (first  term  of  eq.  (B2))  is 
obtained  by  a  Taylor  series  expansion  at  station  n  evaluated  at  the  points 

^n  -  and  ^n  +  ij.  Thus, 
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"\ 


ae\  =  k/de\  H  ±L  Se\  +  (^)2  &L4  &e\  _  <£-4  M) 

J  !  U in  2  SlV  dS/n  8  S|2\  ^/n  48  8|3V  Hjn 


n~2 


>  (B3) 


$±\  +£L  i_4  <£-/ k  ^L(k  |£) 

8  8| 2  \  8£  /n  ^8  ^  V  /  n 


v/M]  =  k(—  ]  +tfLiI_[ki^]  +  H^d_il_[k 

UJ  !  2  ae\  8ijn 

x  'n  +  i 
2 


These  equations  are  then  solved  for  the  second  derivative  and  yield 


i  At  80 


x8|  bin  l 


i/k 


S|4+§ 


i/k  ^  +  ^(k  M' 

24x  Sg3l  ae 


8Mn-I 

2 


(b4) 


n 


or 


1  S  [k  d0 


^n+l  “  -^n 


2  al  ^  a^  jn  ~  1 4n>n+1  *  "  kn~1>n  1  ) 

which  is  correct  to  order  I2*  Note  that  x  A £  =  7. 

The  finite-difference  form  of  the  first  partial  derivative  in  the  second 
term  of  equation  (B2)  is  obtained  by  a  Taylor  series  expansion  about  station  n 
evaluated  at  n  -  1,  n  +  1,  and  n  +  2.  Thus 

e  .  .  8  .  J*£\  +  Mf/afa'l  - 1 + -  .  . 

0n-l  “  en  t*[nl*a^-\glt2)a  6  L j)  *  l,5Vn 


Tn  *  Tn-l/ 


0n+1  -  0n  +  ^(§)n 
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2  W  /n  6  W5/n  WVn 
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This  set  of  equations  is  solved  for  the  first  derivative  and  yields 


I  |f  =  ^(6Tn+l  -  3T„  -  2Tn_l  -  Tn+2)  + 


(Aj)3/'a40 


12 


Therefore,  the  expression 


•^Tn+p  -  3Tn  ~  2Tn-l  ”  ^n+2) 


is  an  approximate  derivative  at  station  n  (2  ^  n  =  i  -  2)  correct  to  terms  of 
order  l3. 

At  station  n  =  i  -  1,  the  series  is  evaluated  at  n  -  2,  n  -  1,  and 
n  +  1  since  there  is  a  discontinuity  at  station  i.  This  procedure  gives 


t.x  ■  U211  +  3Ti-i ' 6Ti-2  +  ti-3) 


which  is  also  correct  to  terms  of  order  I-3. 

Substituting  equations  (B4)  and  (B7)  into  equation  (B2)  gives  the  finite- 
difference  form  of  equation  (B2)  as 


■^n-1  “  ^n 


cn-l,n  1  “  ^n,n+l 


Tn  “  ■*'n+l 


+  |(6Tn+l  "  3Tn  -  2Tn_!  -  Tn+2)  iccp  +  “p5p  +  ~~ “  “c)° 


+  ZF  "  PV  AT 


(2  ^  n  <  i  -  2)  (B9) 


and  by  using  equation  (b8)  instead  of  equation  (B7),  the  equation  is 


k  Tn-1  ~  Tn  _  k 

Kn-1, n  ,  ■K; 


■^n  "  -^n+1 


n, n+1  ~ 


+  7 (2Tn+l  +  3Tn  -  6Tn_p  +  Tn_2^  mcCp  +  mpCp  +  — 


n  -  l/  V 


Tt^7 "  fflc  cp 


+  ZF  =  pcpz  aT 


(n  =  i  -  1)  (BIO) 
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Uncharred  material  ( i  +  1  %  n  1?  i  +  ,j  -  l).-  The  transformed  equation  for 
the  second  layer  (eq.  (24b))  is 


(: 


A  4*  +  JL  V!fP(i  _  nMl  =  p.c*  A1 

c.)2  A )  x1  p’  -  Pu  9  cp  at 


(Bll) 


The  second  partial  derivative  is  found  as  in  the  previous  section  and  is 
given  by  the  following  finite-difference  expression  (correct  to  order  ^2). 


1  _djk»  A 


1 L' 


Vl  -  Tn  ,  . 


t  •  Tn  ~  Tn+l' 

(x-)2^r  hl  =  i'  "  n>n+i  i *  j 


(B12) 


and  as  in  the  previous  section,  the  first  derivative  is  obtained  by  a  Taylor 
series  expansion  about  station  n 


1_  Ml  =  JL 

\x’  &£/n  " 
and 


?(6Tn+i  -  3Tn  -  2Tn_!  -  Tn+2)  (i+l^n^i  +  j-2)  (B13) 


f  1  _  JL 


P^n+l  +  3Tn  -  6Tn_i  +  Tn_2)  (n  =  i  +  j  -  l)  (Bl4) 


These  expressions  are  again  correct  to  order  Z^1. 

This  procedure  gives  the  finite-difference  form  of  equation  (Bll)  as 

rp  _  rp  m  _  rp 

V »  n-1  xn  ,,  I  xn  n-1 

&n-l,  n  ^  i  -  n+1  ^  i 

+  (^■'V  T  'VP  orp  -  rp  ^pP  <""P  I  I J  •  ^Tn 

+  6  V  n+1  3Tn  "  2Tn“1  ■  Tn+2j - J -  -  P  cVl  At " 


(i+l<n<i  +  j-2)  (B15) 
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and  as 


■^n-1  ”  -^n 


n-l^n  1 1 


k' 


Tn  “  •'•n+1 


n) n+1  i ' 


+  i(?rn+1  +  JTn  -  at!  +  M1  - f  —  -  »’p’«p  Is 


(n  =  i  +  j  -  l)  (Bl6) 

Insulation  (i+j  +l^n^i+j  +  m  -  l).-  The  differential  equation  for 
the  third  layer  is  (eq.  (3)) 


8  L«  de"\  _  t.  ..  Be 
5T  ay  j  -  p  °p  at 

The  second  derivative  is  given  by  the  expression  (correct  to  order  l 2) 


(B17) 


l.[k"  Mil  -  k" 

aor\  ay  Jn  '  n-1’n 


Tn-1  ~  Tn 


Tn  “  Tn+1 


n,  n-1  ^  » 


(B18) 


Therefore,  the  difference  equation  is. 


,  _  H'l .  k 


vn-l,n  2" 


Tn  -  Tn+q  _  iti  „  „  ATn 
n, n+1  2»  “  1  p  CP  At 


(B19) 


Boundary  Stations 


First  station  (|  =  0).»  Some  difficulty  is  encountered  in  finding  an 
acceptable  finite-difference  equation  for  the  first  station.  The  first  approach 
that  is  used  to  obtain  such  an  equation  is  relatively  straightforward.  First, 
the  second  partial  derivative  of  equation  (B2)  is  expanded  in  a  Taylor  series 

2 


about  the  point  |  =  0  and  evaluated  at  |  =  —  • 


i_  al 

x2bi{ 


k 


Tc 


Tn 


1,2 


M\ 

^A=o 


(B20) 
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Substituting  equation  (B20)  into  equation  (B2)  gives 


x\*h=o 


+  —  ^mcCp  +  ApCp^  +  F  -  pep  — —  (B2l) 


ATq 

AT 


The  transformed  equation  for  the  boundary  condition  at  the  first  station 
(eq.  (25a))  is 


^•aero  ffel 


_  k  Se 


+  mc 


■where  q, 


aero 


x  d£ 

is  given  by  equation  (25b). 


Si  0  -  Tx  11  Hc  +  Ahc  I  -  Ah 


(B22) 


The  exact  value  of  fi  —  ]  from  equation  (B22)  is  then  substituted 

lx  d§/|  = 0 


into  equation  (B2l)  to  yield 


%2 


T2  ~  T1 


aepTJ  +  mc 


(e  -  Tp^Hc  +  Ahc)  -  Ahc 


-  1. 


•aero/ 


1  - 


:(m, 


2kVmcCP  +  mPCP 


+  I  F 


l 

pCP  2  aT 


(B23) 


The  difficulty  encountered  in  attempting  to  use  this  equation  is  due  to  the 


term  1  -  —  ^mcCp  +  ApcpJ  •  This  term  must  be  approximately  equal  to  1  since  l 

is  an  approximation  to  an  infinitesimal*  However,  since  k  is  typically  of 
the  order  10“^  or  smaller,  l  must  be  taken  so  small  that  the  computer  time 
required  for  an  acceptable  solution  is  unreasonable. 

To  alleviate  this  problem,  a  different  expression  for  the  second  first 
derivative  of  equation  (B2l)  is  used.  This  expression  was  found  by  expanding 
the  temperature  at  station  1  in  a  Taylor  series  evaluated  at  stations  2,  3, 
and  k  which  gives 


i/ae \  = 

xV^A=o  61 


(-llTq  +  18T2  -  9Tj  +  2T| 1 


(B24) 
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Hence  the  finite-difference  equation  for  station  1  becomes 


*1,2  +  ^ero  -  <*1*1  -  %rs(Tl  -  Tl)(Hc  +  ac)  -  ^ 


+  ^(-11T1  +  l8T2  -  9T3  +  2T4|)  (mccp  +  rificpj  +  I  F  =  pcp  |  jji  (B25) 
Note  that  the  exact  expression  for  [i  ^-)  from  equation  (B22)  is  still 

V  ^/|=0 

used  for  the  first  first  derivative. 

Interface  (|  =  1,  C  =  0).-  The  transformed  equation  at  the  interface 
boundary  \  =  1,  £  =0  is  (eq.  (26b)) 


k  de\  •  ak  fit’  de' 


'  *  0 


(B26) 


The  finite-difference  form  of  equation  (B27)  is  found  by  a  simultaneous 
solution  of  equations  (B2)  at  |  =  1,  (Bll)  at  £  =  0,  and  (B26).  In  this  way, 
the  interface  equation  is  satisfied  both  at  the  interface  and  also  immediately 
on  either  side  of  the  interface.  To  do  this,  equation  (B2l)  is  rewritten  for 
|  =1.  Thus, 


k  de\  _  k  xi 

„x  atA-o. " ' 1-1,1  1 


Ti  ~  Ti-1  +  1  (be 


M^/i=i 


mplcp  + 


cpP 

_____ 


7  7 

+  IF-  ^pIaT 


(B27) 


Then  by  using  a  Taylor’s  series  expansion  about  station  i,  evaluated  at  i  -  1, 
i  -  2,  and  i  -  equation  (B28)  can  be  written  as 


k  _  v  i 

-  i-1^i  l 


12V  P  P'  -  P 


(llTi  -  l8Ti_! 


+  9Ti_2  -  2Ti_3)  +  |  F  -  pcp  1 


l  ^i 


(B28) 
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A  similar  procedure  using  equation  (Bll)  at  £  =  0  gives 


!„t 


rk'  de'\  _  ],!  Ti+1  ~  Ti  mp  /  CPP 


§T/£=o  =  "kl’1+1  T'  l5\p^Ty 

\  .  7  t  AT-,* 

-  9Ti+2  +  2Ti+3j  +  P,cp  *2  AT  (B29) 

Finally,  substitution  of  equations  (B28)  and  (B29)  into  equation  (B26)  gives 


-iiTi  +  i8ri+1 


i  Ti+1  “  Ti  T-l  -  Ti_i  .  I /_  Pcp 

ki,i+l  ^  7 . .  +  + 


_1_ 

12 


(llTi  - 


i8t 


i-l 


+  9Tq_2  -  2Ti_3) 


\ 

]  +  Cpp' 

r~ 

1 

V 

p'  -  p 

12 

;(-llTi  +  l8Ti 


i+1 


9Ti+2  +  2Ti+3^ 


-  AhT 


+  iF  =  foe  1  +  p«c*  ±A— 

2  V  p  2  P  P  2  /At 


(B30) 


Back  surface  (2  layers,  t  =  1).-  For  a  two- layer  system,  the  transformed 
b oundary  equation  is  (eq.  (27)) 


kVae*  _  _ 

x'U  ^=1  "1+J 


C^a.^  ~  +  S 


(e '  -  Ti+J)Awf  ^f 


+  aei+j 


(e'f  -  4 


(B31) 


The  first  term  of  equation  (BJl)  is  determined  by  evaluating  equation  (Bll) 
at  £  =  lj  thus, 


fk1  de! 


Ti+j  ~  Tj+j-l  ,  l '  m'i+j 


”ki+j-l, i+j  2’ 


,  m 


\X’  *  ^=1 

Substituting  equation  (B32)  into  equation  (B3l)  gives 


P  c. 


P  2  At 


(B32) 


ki+j-l,i+j  ji 


^  ^  -  «i+J(«4+J  -  4) 


1  (Ti+j  -  Ti+j^A^i  Ahp  (p  cp  2  + 


^i+j 

At 


(B33) 
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Interface  (3  layers,  £  =  l).-  For  a  three-layer  system,  the  transformed 
b oundary  equation  at  £  =  1 is  (e q .  ( 28b ) ) 


'k 1  1  \  P>  .  _ 

vx'  k=i  i+J 


Se '  _  ( Sb 


Sy 


(B?4) 


'y=x0+x0 


The  first  term  of  equation  (B34)  is  given  by  equation  (B32).  The  last 
term  of  equation  (B34)  is  evaluated  by  expanding  the  second  derivative  of  equa¬ 
tion  (Bl8) .  Thus, 


Lr*S-\  -  k" 

V  ^40+x-"  i+^i+J+1 


Ti+j+l  "  Ti+j  „  «  l"  ^i+j 


+  P  °p 


2  At 


(B35) 


Substituting  equations  (B32)  and  (B35)  into  equation  (B34)  gives 


1:i +.]-!,  i+j 


Ti+j-l  -  Ti+j  ■■  T<-  -  T< 


-  k. 


i+j, i+j+1 


-i+j  ~  i+j+1 


f  7  • 

p'c'  —  +  p"c"  h~  +  c 


vAT., 


_  i+j 

p  ~2  1  “i+j  I  At 


(B36) 


Back  surface  (3  layers,  n  =  i  +  j  +  m).-  For  a  three-layer  system,  the 
back-surface  boundary  condition  is  (eq.  (29) ) 


-k"  =  Ci+j+m  ^  +  s(e"  -  Ti+J+m)A/f  Ahf  +  ue 


Sy 


i+j+m 


o")4  -  a* 


(B37) 


Equation  (B37)  is  solved  simultaneously  with  equation  (B17),  with  the 
previously  used  approximation  to  the  second  derivative,  to  yield. 


v" 

i+j+m-1, i+j+m 


rn  _  rn 

i+j+m-1  i+j+m 


^i+j+ml^i+j+m  “  ^b)  +  ®(Ti+j+m  “  ^i+j+m^^f  ^f  +  (p"cp  2  +  ^i+j+mj 


\AT 


i+j+m 

At 

(B38) 
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SIMPLIFICATIONS  TO  REDUCE  COMPUTING  TIME 


Computing  the  Time  Interval 


The  primary  difficulty  encountered  in  using  the  program  is  the  time 
required  to  obtain  solutions  on  the  digital  computer.  The  equations  are  inte¬ 
grated  stepwise  over  time,  and  the  maximum  time  interval  for  which  the  solution 
is  stable  with  the  explicit  formulation  used  is  (ref.  25), 


At 


pcp22 

~2k 


where  the  properties  of  the  layer  giving  the  minimum  At  must  be  used. 

Theoretically,  the  initial  value  of  l  may  be  zero  in  the  physical  sense; 
however,  a  finite  initial  char  thickness  is  required  for  the  calculating  pro¬ 
cedure  used  in  this  program.  The  initial  value  of  the  char  thickness  which  is 
usually  chosen  is  extremely  small.  As  the  thickness  of  the  char  layer  increases 
due  to  the  pyrolysis  at  the  interface,  the  preceding  stability  equation  allows 
the  digital  computer  to  use  a  larger  time  interval  in  its  stepwise  calculations. 
Likewise,  when  the  second  material  pyrolizes  to  such  an  extent  that  its  l 
becomes  very  small,  the  preceding  equation  now  decreases  the  time  interval  so 
that  the  finite-difference  equations  remain  stable.  This  time  interval  is  only 
allowed  to  change  whenever  the  digital  computer  can  recognize  this  change  and 
continue  to  print  out  answers  at  the  specified  print-out  times.  Also  the  com¬ 
puter  recognizes  specified  minimum  values  for  both  x  and  xT  so  that  com¬ 
puter  time  is  not  wasted  when  either  x  or  xT  approaches  zero. 


Options  for  Thin-Layer  Configurations 


In  some  cases,  it  is  necessary  to  run  calculations  where  the  stability 
time  interval  At  becomes  very  small.  Therefore,  to  decrease  computing  time 
by  at  least  one-half,  the  computer  recognizes  input  quantities  b  and  b 1 . 
When  x  1*  b  and/or  x1  ^  bT,  the  interior  stations  are  no  longer  considered 
and  the  entire  layer  becomes  two  half  elements  of  thickness  x/2  and/or  x!/2- 
In  this  procedure  the  following  finite-difference  equations  are  utilized. 


When  x  %  b,  the  equations  at  the  first  station  (eq.  (B25))  become 


2kl,i 


Ti  "  T1 


+  2  kero  "  CTelTl  “  ®c 


s  T-l  -  TJIHc  +  Ahc 


Ah, 


+  (Ti  “  TlX“ccp  +  “pCp)  +  xF  =  pepX 


ATj 

At" 


(Cl) 
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The  equations  for  the  interior  stations  in  the  char  layer  are  no  longer  con¬ 
sidered  when  x  ^  b.  The  equation  at  the  interface  (eqs.  (BjJO))  is  as  follows: 

When  x  b  and  x '  >  b  ' 


2k 


i,i+l 


^i+1  “  ^i 

r 


2k 


M 


+  xF  + 


-  T 


x) 


+ 


*  . 

CpP* 


P'  -  P 


(-llli  +  18T1+1  -  9Ti+2  +  2T1+3)  •-  2Ahpl  =  (pcpx  +  P'cpZ 


When  x  ^  b  and  x '  ^  b ' 


,  Ti+i  -  T±  Ti  “  T1 

— P - a!,i  — —  +  **  +  "I 


pc. 


c^  + 


p 


p  p'  -  p 


-)(' 


cpp  , 

+  - (T 

P 


?ry(Ti+,j  "  Ti)  -  *  (<*V  +  p,cPx')at 


A^l 


')^  <C2> 


(03) 


When  x  >  b  and  x  *  ^  b ' 


2k 


T 


i,i+j 


i+j 


T, 


-  2k. 


x  1 


1-1,  i 


i(iiTi  -  1%.! 


+  9Ti_2 


+ 


CpP' 

P’  -  P 


(c4) 


When  x'  <  b',  the  interior  stations  in  the  second  layer  are  not  con¬ 
sidered.  The  equations  at  the  back  of  the  second  layer  (eq.  (B33))  for  two 
layers,  (eq.  (B38)  for  three  layers)  are  as  follows: 


When  x  ^  b’,  equation  (B33)  becomes 


T,  -  T 


i+j 


x 1 


ae 


i+j 


s(th-3  -  Ti+j)*f 


AT 

P'CP  T  *  Ci+j)-#4  <«> 
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and  equation  (B38)  becomes 


ki,i+j 


Ti  ~  Ti+.1 


"  ki+j^i+j+l  z 


Ti+,j  “  Ti+j+l 


=  fp'cp  ^  +  p"c;  ^  +  ci+j )-£T- 


(C  6) 


The  computer  uses  the  standard  equations  whenever  x  >  b,  and/or  x'  >  b'. 
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